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Abstract. We use the cavity method to study parallel dynamics of disordered Ising 
models on a graph. In particular, we derive a set of recursive equations in single site 
probabilities of paths propagating along the edges of the graph. These equations are 
analogous to the cavity equations for equilibrium models and are exact on a tree. 

On graphs with exclusively directed edges we find an exact expression for the 
stationary distribution of the spins. We present the phase diagrams for an Ising model 
on an asymmetric Bethe lattice and for a neural network with Hcbbian interactions on 
' ^ I an asymmetric scale-free graph. 

For graphs with a nonzero fraction of symmetric edges the equations can be solved 
, for a finite number of time steps. Theoretical predictions are confirmed by simulation 

results. 

Using a heuristic method, the cavity equations are extended to a set of equations 
^ I that determine the marginals of the stationary distribution of Ising models on graphs 

with a nonzero fraction of symmetric edges. The results of this method are discussed 
, and compared with simulations. 
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1. Introduction 



^ I Many problems in different research fields are based upon tlie interaction of units 
through some underlying graph. Some examples are: gene expressions in boolean 
networks [1], agents competing for some limited resources [21 [3], interactions between the 
decoding variables in low-density parity check codes [1], interactions between humans 
on a social network [5], the analysis of phase transitions of a spin glass [6]. 

To calculate statistical quantities on a given graph instance, one can use the 
cavity method [7]. This method is based on the assumption that for sparse graphs 
the neighbouring spins only depend on each other through their direct interactions. 
A similar method, known as the sum-product algorithm, is used in information theory 
and artificial intelligence, see [8] for a tutorial paper. Examples of problems investigated 
with the cavity method are: the characterisation of the set of solutions of optimisation 
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problems on random graphs [1], the calculation of the eigenvalue spectrum of sparse 
random matrices [10] and the solution of the minimum weight Steiner tree problem [11] . 

For many problems the stationary distribution of the spins is not known, e.g., neural 
networks with asymmetric couplings [12], the minority game [2], ... One would like to 
generalize the cavity method to treat the internal dynamics of these models. One of 
the standard approaches here is the generating functional analysis first discussed in |13j . 
The stationary solution of the minority game was found using this method [2] . In work in 
progress [H] the parallel dynamics on graphs is studied using the generating functional 
method of [15]. The same authors used this method to analyze the evolution of a 
decoding algorithm on a sparse graph [16]. Another successful method is the dynamical 
replica analysis which was applied on finitely connected systems [T5| [T7] . Recently, 
sequential dynamics of an Ising spin glass on a Bethe lattice with binary couplings was 
solved using a cavity- like approach [T8] . 

In this work we apply the cavity method to solve the parallel dynamics of models 
on random graphs. The aim thereby is twofold. First, we want to generalize the effective 
dynamical equations that solve the dynamics on a Poissonian graph, given in [15], to 
random graphs with a given degree distribution. This generalisation is important when 
we want to solve the dynamics of models on a given graph structure. We use these 
equations to solve the dynamics and to find the stationary distribution of an Ising 
model and a neural network model, both on a fully asymmetric graph. We discuss how 
the correlations in the indegrees and outdegrees influence the performance of the neural 
network. The second purpose of the paper is to extend the cavity equations to graphs 
with both symmetric and asymmetric couplings. For this we need to find the stationary 
state of the dynamics. This is possible when we neglect the correlations in time of the 
stationary distribution. We discuss how good this approach can predict macroscopic 
observables of Ising models with bond disorder or with fluctuating connectivities. 

This paper is organized as follows: In section [2] we define the necessary quantities 
and derive the effective dynamical equations for the single site marginals on a given 
graph instance. We take the average of these over a graph ensemble in section [3l In 
section H] we specify the dynamics of the spin models. We derive the equations for 
the distributions of single site marginals in section [H The evolution of macroscopic 
observables obtained from the theory is compared with simulations in section El We 
discuss the phase diagrams for an Ising model on an asymmetric Bethe lattice and for 
a neural network with Hebbian interactions on a scale-free graph in sections [7] and M 
respectively. In section [9] we derive an algorithm that calculates the single site marginals 
of the stationary distribution on graphs with arbitrary symmetry. A discussion is given 
in section [TOl 
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2. Dynamics on a given graph instance 

2.1. Some Definitions and Notations 

We consider models defined on a given graph instance G = {V,E), with V and E 
respectively the set of vertices (or sites) and the set of edges. We limit ourselves 
to simple directed graphs G determined by a connectivity matrix C, with elements 
[C]-j = Cij e {0, 1}. When = 1 and cji = the graph has a directed edge from the 
z-th site to the j-th site. When Cij = Cji = 1 there is an undirected edge between i and j 
and when Cij = Cji = there are no edges between them. We define the sets E, E'^ and 
^sym through: E = {{ij) eV x V\cij = 1}, E'^ = eV x V\cij = l,Cji = 0} and 

^sym ^ {{i,j) ^ V X V\cij = 1, Cji = 1}. We study the evolution of Ising like models of 
n-replicated Ising variables o"j(t) G { — 1, 1}", with i = 1, . . . N and t the corresponding 
discrete time step. The dynamics in discrete time is defined by a transition probability 
W{cr{s)\(T{s — 1)), from the state cr(s — 1) = {<Ji{s — 1), a2{s — 1), ■ ■ ■ , o"Ar(s — 1)) on 
the (s — l)-th time step to the state cr(s) on the s-th time step. We consider transition 
probabilities W of the form: 

N N 

W [cTis)\cTis -iy,0] = l[w [a.(s)|(T(s - 1); 6.(3)] = l[W [a.{s)Ms)] . (1) 

i=l i=l 

The n-dimensional local field hi{s) is defined through 

h,{s) = J2 hj^^Ms ~ 1)) + 9,is) , (2) 

where the field hj^i{aj{s — 1)) quantifies the influence of the spin on site j on the 
spin on site i and 6i{s) is an external field. We used 9^'" for the neighbourhood of 
all the vertices that influence i directly, i.e. (9™ = {j E V\cji = 1}. We will also use: 
^out ^ ^Y\c.. = i^^ g.^ gin y ^out ^sym ^ gin p ^out_ rpj^^ probability to have 

the path cr*"--* = (cr(to), ■ ■ ■ , cr(^)), from time step to to time step t, is given by 

Pt...t (a*-t|0*o+i..*) =( II W [ct{s)\ct{s - 1); e{s)]\ {a{to)) , (3) 

\s=to + l J 

with Pt^ (cr(to)) the probability distribution of the spins at time step to. 



2.2. Dynamical Version of the Cavity Equations 

Using the cavity method, see |^, it is possible to solve the parallel dynamics on 
graphs. The cavity graph G^*-* is the subgraph of G where the i-th vertex and all of 
the interactions with its neighbours are removed. We write the following relationship 
between a path probability Po..t on the graph G and the probability Pq \ on its related 
cavity graph G*^*^ 

Po..t K'l^^-*) = Pill {a^-'\e'-' + C«'^-*) ^ [^^(^)l'^(^ - ^^(^)]) Poi^^m ■ 

(4) 
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In equation (j4j) we introduced an extra field (j^\s), representing the influence of the i-th 
spin on its neighbours j G d°^^: 

Cf{s) = e,,h,^,{a,{s-1)). (5) 

The prefactor etj determines whether the edge is symmetric or not: etj = 1 for undirected 
edges and e^j = for directed edges. We took a factorised initial distribution Pq: 
Po(c(0)) = Hill Po(o'j(0)). The single site marginal Pi,o..t is obtained by summing Po..t 
in (jl]) over all paths cr° '* with j ^ i ■ In general, we will use the notations 

0"5 = (crn,0-j25 ■ ■ ■ ,0-i,s|) > (6) 

Ps{ar\e'-') = ^('^"■lO, (7) 

with S a set of indices: S = ■ ■ ■ ,i\s\} , where \S\ denotes the size of the set S. 
Within this notation Pa-.o..t is the joint probability of the paths on the neighbours of i. 
When we sum over all paths crf'^, with j ^ i, on the left hand and right hand side of 
(HD, we get 



r0..t s=l 
9i 



In the sequel we drop the subscript i in the argument of Pi.o..t- Now we make the 
Bethe-Peierls approximation: i.e. we assume that the spins in the neighbourhood di of 
i become independent when we remove the i-th spin: 

pS, (j!^ ■ • • . i«" + «"'■'■■') = n ^S... ('■° 'I*]-' + ^«c!""') . (9) 

with di = {jl, ■ ■ ■ ,j\di\}- In ([9]) we took 9j = when j ^ di U i. We substitute Qj in 
(E]) to get the following set of recursive equations 



0..t-l \j(zSin\i 
/ t \ 



l[w[a{s)\hf\s)]po{am] , (10) 

for the path probability Pf^ j on the graph G*^^-*, with i G di. To derive ( |TOl) we used 
-^i'a.t ~ -^ia.f '^'^^ I -E I -equations ffTOl) determines the -probability distributions 

p/q^ J ((T*^"*|^j^ '*) at time step t as a function of the [E'l -probability distributions 
-Pi''a.i-i('^°"*~^l^i^ *^^) previous time step t — \. In equation (fTOl) we only need 

to take the product over j G (9™ because the flelds hf\s) depend only on (J QiYi . We 
call the equations ffTOl) the dynamical cavity equations analogous to the static equations 
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(IB.14P . The main difference is that (11 01) are recursive equations of probabilities of paths 
propagating along the graph, while in fIB.Mp messages are propagated that determine 
the marginal probabilities of the stationary distribution. Just like the static equations, 
the set of equations (fTOj) is exact on a tree. To find the marginal distributions Pifi,,t on 
the original graph G from the cavity distributions, we need to combine equations (IHl) 
and dlD: 

^O..t-l \j Qin 

i 

\{W[a{s)\K{s)]p,{a{Q))\ . (11) 

Equations ( fTTi) are the dynamical versions of the set of equations (IB. 151) . The initial 
problem of finding the single site marginals Pj,o..t from the A^-site probability Pq..* has 
a computational complexity 0{2^). The set of equations ( fTOl) and ( fTTl) has a linear 
complexity 0{N) in the system size and an exponential complexity O (2*) in time which 
makes the dynamics solvable for a finite number of time steps. 

The cavity equations simplify a lot when the graph is fully asymmetric. In this 
case we can set e^j = in equation (ITOl) . Therefore, the equations only have to be 
solved for 6^''^ = 0°"*, where 0°"* is the null vector. Moreover, because eij = the 
self-coupling disappears in ffTOj) . We can thus sum on the left and right hand side of 
(HO]) over (or,(0), a,(l), ■■■,a^{t~l)) to get 



Equation ( |T2l) describes a Markovian dynamics. 



3. The Ensemble Averaged Distribution of Paths 

We calculate the average of equation ffTOl) over all links in the graph, i.e. all a & E. The 
graph is drawn from an ensemble of graphs Q. We look at ensembles where the typical 
graphs have a local tree structure and the degrees on different sites are uncorrelated. 
An example is the Poissonian ensemble Qp defined in ( lA.ip of Appendix A The degree 
distribution is defined through a histogram as 



In equation (fT3l) we use the following notations: the indegree kf" = \df^\, the outdegree 
^out _ I ^out I ^YiQ symmetric degree kf"^ = \ 9™ fl df^^ \ . For N oo the dynamics 
of Ising models on typical graphs drawn from such ensembles depends on the degree 
distribution (fT3l) . We define P as the average of the path probabilities P- f over all 
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directed edges of E: 



p"/ 0..A _ ^(i/)eEd i,o..t 



(14) 



— sym (£) 

The average probability mass function P is defined as the average of P// over all 



links belonging to an undirected edge 

pW (^O..t\nl..t\ 



I E'sym I 



(15) 



When we use the property that the spins in the neighbourhood of i are uncorrelated, 
we can write 



^0..*-l|^«,l..t~l 



(16) 



It is useful to focus on a specific example. We consider fields of the type hj^i{s — 1) = 
Jji(7i{s — 1), where the interactions strengths Jij are i.i.d.r.v. drawn from a distribution 
R{J). When we take the average of the update equations ffTOl) according to the 
definitions (fT4l) and (fTSi) . and use ( fT6|) we find the recursive equations for the averaged 
probability mass function of paths. These recursive equations are given by: 

oo oo min{A;°"t,A:'") 



Cout C, 



sym 



i-1 



po(a(o))n^ 



s>0 



a(s + l)|e(s)+ J] Ji'ae{s) 



0<^'<fc"> 



:i7) 



and 



-sym 



£=1 ,tO-'-1 



t-1 



s>0 



a{s + l)\9{s + 1) + E -''^'^^'(s) 



0<^'<fc'n-l 



We introduced the average connectivities Cgym = Xlfci" fco^^t ^sym p(/i;^°, Z;;'^-^™)/;;^'^™ and 

— TiTi-in tout ii,sym\iLout 

'-'Out / jf^in f^out fcsym pyrh ^ , J fb 



■real , 



The averaged probability mass function P (cr ) over the marginals Pi{a 
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defined through F = J2^P^i(^°^^')/N, can be calculated from ffTTl) : 



t-l 

poKo))n^ 



s>0 



(19) 



a(s + l)|e(s + 1) + Y Je'^^i'i^) 

The Markovian dynamics of spins defined in ([T]) is thus reduced to an effective non- 
Markovian dynamics of one single spin given by the recursive equations f|T7|) . f|T8|) and 
f|T9|) . Equations analogous to f|T7|l and f|T8|) were derived in [16] in the context of LDGM 
channel coding using the generating functional analysis. 

For fully asymmetric graphs, see (fT2|) . we remark that Pj^f{a) = Pi^t{c), but the 



averages, (a) = P'>^^{a) and p^""" (o") = Pi t{a), over, respectively, the links and the 
sites are different. Indeed: 



irrcal 



5^Kn^ n f dJeR{Je)YPU^^) 

t 

Pi^mU^Hsms)] , (20) 
rr\^) = J2p(k"') n /rf^.i?(^.)E^'-iM^'(^(o))ri^[^(^)i^(^)] ' (21) 

with c(/e"^) = Efcoutp(/e°"N^'")^°"* and Cout = Efcout p(A;°'^')A;°"^ 
4. Examples of Dynamics 

In this section we define the type of dynamics we study by specifying the form of the 
transition probabilities M/^ [(t|/i] used in equation ([1]). 

4.I. Glauber Dynamics 

We consider Glauber dynamics for an Ising model with n = 1, i.e. ai G {—1, 1}. Every 
spin (Tj(t) evolves under the influence of the field hi{t — 1) with a transition probability 
Wg{ai(t)\hi(t)) defined through: 

The parameter /3 is the inverse of the temperature T. It is possible to implement the 
dynamics defined by (l22l) and ([H) with the heat-bath algorithm [19] . When the graph is 
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fully symmetric detailed balance is satisfied and the Hamiltonian is given by equation 

As an example we apply the recursive equations ( |T71) . ( |T8|) and ( |T9l) to a Poissonian 
ensemble Qp defined by the probability Pp{C; c, e) of the connectivity matrix C: 

Pp(C; c, 6) = n {^Sic^r, 1) + (l - ^) 0)) 

i<j 

n c,.) + (1 - 6) (^5(c,,; 1) + (l - ^) 5(c,,; 0))) . (23) 

The parameter c controls the number of edges while e controls the fraction of symmetric 
edges in the graph. The S{a; b) are Kronecker delta functions. For N — > oo, the typical 
degree distribution of a graph drawn from this ensemble, Pp(A;'°, A;°"*, k^^"^), is given by: 

Pp(A;^", F^-) = exp [-C6] exp [- (1 - e) c] ^ 



^sym] M (A;'° - fc'^y'")! 

For the derivation of (l24l) we refer to Appendix A[ Substitution of (12^ in (fT8|) and (fT9l) 
leads to, 

k 



fc>o ■ \f=l 
II [eP''"K-*-V^^°-*"') + (l-e)F""K-*"'|0°-*-' 



£=1 ^O-.t-l 



^^^^'^^ l[ 2 cosh [/5 (^(s + 1) + Eo<.<. ^.'^H^))] ■ ^ ^ 



For the Poissonian ensemble equation (]A.9P is valid such that P (cr° ■*) = 
p^y™^^o..t|gi..t-i^ _ p^°^ (cr*^- *). We find equation (l25l) which is identical to the main 
result of reference |20] derived by calculating the generating function. Hence, we 
conclude that the equations f|T71) . f|T8l) and f|T9|) are consistent with the results found in 

Eol. 



4-2. Coupled Dynamics 

We define a dynamics of two sets of N spins under the influence of the same thermal 
noise. We thus have n = 2, i.e. the dynamic variables are {ai(t) , Ti(t)) G {—1,1} x 
{— 1, 1}. The spins (crj(t), Ti(t)) feel only the influence of their neighbouring spins (aj, Tj), 
with i G 5™, through, respectively, the flelds hi(t) and gi(t). The flelds /ij(t) and gi{t) 
depend, respectively, only on the sets cr(t — 1) and T(t — 1). The spins {ai{t) , Ti{t)) 
evolve according to [{ai{t),Ti{t))] hi{t), gi{t)] : 

Wc [{a,T)\h,g] = 6 {a; -r) \rh - rg\ 
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-6 (a; r) (1 - \rh - rg\)Q{rh - Tg) 



+5 (cr; r) (1 - \rh - rg\)<c){rg - th) 



5{a-l)- 



+ 5{a- 



where 6 is the Heaviside step function and the weights and Vg are given by 

exp {i5K) exp {jig) 



rh 



(26) 



(27) 



2 cosh (ph) ' ^ 2 cosh {/Sg) ' 

Equation fl26l) can be simulated using a heat-bath algorithm where at each time step we 
choose the same random numbers for both set of spins cr and r. A more compact form 
of W, is: 

Wc[{a,T)\h,g] = 6 {a;-T)\rh- rg\ 

+ 6 (a; t) e(r, - r,) [5(a; 1) r, + 5(a; -1) (1 - r,)] 

+ 6 {a; r) 0(r, - r,) [6{a; 1) r, + 5(a; -1) (1 - r,)] . (28) 

When the thermal average of the distance between the paths <T(t) and T(t) does 
not converge to zero for t —>■ oo, even when the initial distance between cr(0) and t(0) 
is very small, the system is in a chaotic phase. We use the transition probability Wc 
to determine the phase transitions to this chaotic phase. Chaotic behaviour has been 
studied in [21] for spin glasses and in [22] for neural networks. The coupled dynamics 
(!28|l can not satisfy detailed balance. 



5. The Path Entropy and the Distribution of the Probability Distributions 
of Paths 

The fluctuations of the path probabilities Pj^^^j ^ over all links are given by the distribution 
of the probabilities of the paths which we will call V. They determine quantities like the 
average path entropy S{t). On the basis of the recursive equations for the distributions 
V we discuss in section |9] the stationary solutions of the dynamics. 
The average path entropy is deflned as 



s{t)^-j2Po..t{^n iog(Po..t(^°-*)) 



(29) 



where the bar denotes the average over the quenched variables. With the cavity method 
[7|, we can write 

s{t) = Y^Asnt) -- J2 - E ^'^""'(^) • (30) 

i=l aG^sym ae£;d 

The quantity Ai5f*°(t) is the increment in the entropy S when the i-th site is added to 
the graph G^^^: 

t 
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log (|po,o..i(4-*) [^(^)l^(^)]Po(a(0))j , (31) 

with 

The quantity AS^^'^it) is minus the entropy difference when remove the linlc a from the 
graph G 

^o..t ^o..t 

logfP,o..(^°-*)PS..,(r°-V.,a°-*-i)) . (33) 



The summation over the sites in equation (!30|) can be done when we know the 
distributions of the probabilities of paths V on the graph. 
We define the following distributions 



-P'iP) "-J^^ (34) 

0..t\nl..t\ tdW f„0..t\nl..f' 



^""^^ (P) - — ^ '-JE^\ ' (35) 



■preal ^p)^ _ A^t=l 1 lo-"- ' " V V / // (36) 

When the variables Jij are i.d.d.r.v. satisfies the recursive equation: 



oo 



oo min(fc°"Sfc'") 



fcout>0fcin>0 fcBym^o ^^^^ 



n /" ciJ,i?(J,) /" dPeV'^^iPe) n / ^^^^('^^) / dP.V^'^'iP,) 



with 



0..t-2\ 



o..t-i ... J..t-i £=1 

t-1 



The distribution along symmetric links is given by P^y™ 

psym (p) = ^ 1^ p(fcm^fcsym)fesym 



(37) 



V r>o 2cosh [/5Eo<f'<fcin^£'f^^'(s)i ' 



The Cavity Approach to Parallel Dynamics of Ising Spins on a Graph 1 1 



Yl I dJ,R{J,) I dPiP'^^ {Pe) II I dJeR{Je) I dP^V'^' (P,) 
with 



(39) 



0..t-l .^^ 0..t-l f=l 
1 ' ' fein-1 



11 ) Po(f^(o))ii — - — 1 \nn( I n I -J — 



(40) 



The distribution of the single site marginals Pi on the original graph is given by 
V{P)=Y, Yl Pik'^k^n n / dJeR{Ji) / dP,r'^^{P,) 

^in 

n /" rfJ,P(J,) /" dPeV" (Pe) llsfpia'-') - r^^' {Pe, M e=i..k^^^) ) (41) 
with 

fcsym 

^-^^ {P„ J4..i....n) = E n ^.(^-^-V.^.""-^) 

0..t-l „0..t-l £=1 
1 ' ' fcin 

,=^■^+1 V 2cosh[/?Eo<,,<,i„ J,,a,,(s)J y 



In section |9] we use the equations fl37j) . (!39|) and (14T]) to derive the stationary limit from 
the dynamics. When we compare the equations (13 7p and 0391) with the density evolution 
equations flB.12p we see a couple of differences. Since the graph is directed we have now 
two distributions: one for the probabilities propagating along symmetric edges and one 
for the probabilities propagating along directed edges. Since the equations (137|) and (!39ll 
describe the dynamics of the model they are recursive equations. The computational 
complexity of 0371) and 0391) scales exponentially in time. Equation 04 ip is the dynamical 
equivalent of OB. 131) . 



6. Comparison with Simulations 

In this section we compare the magnetisation m{t) = "^^o.-t o'{t)P'^'^'^^ {a^"^\0^"^), 
predicted by equations f|T7j) . f|T8|) and f|T9|) . with results from simulations. It is difficult 
to develop an Eisfeller-Opper scheme [23] for these equations because the probability 
distributions of paths p^-^™ depend on the fields 6'^"*, such that it is necessary to solve 



The Cavity Approach to Parallel Dynamics of Ising Spins on a Graph 12 



— T/(JoC)= 1.5 




) I , I I I I I , L 

2 4 6 8 

time 



Figure 1. The magnetisation m as a function of discrete time for an Ising model on 
a symmetric Bethe lattice with connectivity C — 3. The interactions are drawn from 
the bimodal distribution (j43|) with p = 0.25. The exact enumeration of the recursive 
equations ([T8|) and (fT9|) (lines) are compared with Monte Carlo simulations (markers). 
The red line is calculated at the critical temperature w 1.3459. 




time 



Figure 2. The time evolution of the magnetisation m at T = 1.8 for an Ising model 
on a Bethe lattice without bond disorder for different levels of asymmetry. Results are 
shown for graphs with fixed indegree fc™ = 3 and a given fixed outdegree = C. The 
exact enumeration of the recursive equations (1171) . (fT8|) and (fT9|) (lines) are compared 
with Monte Carlo simulations (markers). 



f fTSj) for all 2* possible values of 9^ '^. We calculate the first time steps through an exact 
enumeration of the equations (fT71) . f|T8|) and f|T9|) . 

In figures [H and [2] we compare the magnetisation in the first time steps with the 
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results obtained through Monte Carlo simulations. In figure [T] we present the results 
for an Ising model with bond disorder on a symmetric Bethe lattice. The interactions 
J are drawn from the bimodal distribution R: 

R{J) = (^) KJ - Jo) + (^) S{J + Jo) , (43) 

with p the bias in the couplings. In figure [2] we show results for the Ising model on a 
Bethe lattice without bond disorder. Both, the exact enumeration and the simulation 
give consistent results confirming the correctness of equations (ITTl) , ffTS]) and (HM . 



7. The Ising Model on a Fully Asymmetric Bethe Lattice 

The dynamical cavity equations (fTOj) simplify to (fT2l) for fully asymmetric graphs. 
We illustrate this for the Ising model on a graph with a degree distribution 
PB(fc''^,A;°"\ fc'^y'") = 5{k''' - C)p{k'"'^)6{k''y'^). We call the graphs drawn from this 
ensemble asymmetric Bethe lattices. Because now (a) = P^{a), we only need 
to solve the recursive equation (120|) . We discuss this model for two typical dynamics. 



7.1. Glauber Dynamics 



In this subsection we let the spins evolve through Glauber dynamics defined in 
section 14.11 The probability of o"j(s) is given by equation (12^ with a field hi{s) = 



7.1.1. Iterative Equations 

We derive the equations that give the evolution over time of the following macroscopic 
observables: the average magnetisation, the correlation function and the distribution of 
magnetisations. 

We define the magnetisation m{t) through the relation Pt(o") = {m{t)a + l)/2. 
From equation ( l20l) we then get for the evolution of the magnetisation under Glauber 
dynamics 



m{t + 1] 



n 

o<e<c 



EC- 



a limit) 



^ tanh 



(3[e{t)+ 



o<e<c 



(44) 



Jl,J2,---,Jc 



To find the correlation function C{t, t') between spins at time step t and t', we sum 
equation (fT7|) over all spins except a{t) and cr{t'). We get the recursive equation for the 
two time marginal Pt,t'{o'{t),a{t')). When we define C(t,t') through 

Pt,t' {(T, r) = ^ [1 + m{t)a + rm(t') + C{t, t>r] , (45) 
we obtain the recursive equation for the correlation function: 

[1 + a(>m{t) + Tem{t') + aereC{t, t'] 



C{t+l,t' + l] 



n E 

o<e<c Kan 
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^ tanh 



. V o<e<c J 



)<£<C / 



(46) 



tanh \P {9{t') 

o<e<c J } ' Ji^J2,- -,Jc 

When we calculate the distribution of the marginals P/^"* from equation fll2p . or 
equivalently the distribution Wt{m) of the corresponding magnetisations, the following 
recursive equation appears 

c ^ c 



Wt{m) = I \[{dJ,R{J,)) I \[{dm,Wt-i{m,)) 

1=1 1=1 



0<f'<C 



(47) 



The time evolution determined by the equations (H4l) . ( I46l) and ( H71 is confirmed by 
numerical simulations. 



7.1.2. A Stationary Solution 
Using the above equations (14^ . (146|) and (H 
consider the stationary solution m{t + 1) 
shows that m is a solution of 



we can find the stationary solutions. We 
m{t) = m. Substitution of this ansatz in 



m 



n 

0<f<C 



E 



^ tanh 



(4^ 



(3 I ^ Ji>ai 

0<i'<C 

The model has a phase transition between a ferromagnetic phase (F-phase) with 
m > at low temperatures and a paramagnetic phase (P-phase) with m = at high 
temperatures. Because this transition is continuous it is possible to determine the P to 
F phase transition line with an expansion of the right hand side of (148|) around m = 0. 
The critical inverse temperature (3* between the P-phase and F-phase is the solution of: 

l = p2~^^( '"^ I |2r- A;|tanh(/3V|2r- A;|) . (49) 




Equation fjl9|) holds for the bimodal distributions R{J) of the form (j^3l) . Using the 
stationary ansatz q = C(t,t') = C(t — n,t' — n) in ( 146!) we can try to find a phase 
transition between a paramagnetic phase with q = and a spin glass phase (SG-phase) 
with g > 0. Analogously to [20] we find that q = for all temperatures and biases 
p. In figure [3] we show the P to F transitions (solid lines) for different values of the 
connectivity C as a function of the temperature T and the bias p in the couplings . 



7.2. The Chaotic Phase 

Although the spin model studied in this section has no SG-phase, it has a chaotic phase 
(CH-phase). In order to find this phase it is necessary to consider the dynamics of two 
set of spins cr and r that interact on the same graph with the same thermal noise with 
a slightly different initial configuration. The transition probability of the spins (cij, r^) 
is given by ([28]) with hi{s) = Xljeaj" 'Jji'^ji^ - 1) 9i{s) = Zlj-gsi" Jji'^ji^ - 1). 
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Figure 3. The P to F transition lines (solid) for the Ising model on an asymmetric 
Bethe graph are presented as a function of the rescaled temperature T/{JqC) and 
the bias p in the couplings (see equation (|43p ). Phase transition lines for different 
connectivities C are shown. The dashed lines enclose the regions where the model 
behaves chaotically. 



As done in ([7]) it is possible to derive the recursive equations for the single time 
marginals P^{o-,r). We define the magnetisation nit and the thermal average of the 
Hamming distance dt between the sets cr and r through 

1 

- ^ (50) 



4 



, 1 + nita + rritT + (1 - 2dt)aTj 



In the case of a bimodal distribution R{J) we get for dt the recursive equation 

c~i 

dt = 





2 cosh (/5(x+ \y\)) 

1 +pmj_i/(l - dt^i] 



exp (/9(x — 


m 


2 cosh {/3{x - 


-\y\)) 



pmt-i/{l - dt^i) 



5 {x — 2v + n) 



2"-'='"(5 {y-2w + k"' - n) 



(51) 



The time evolution of the magnetisation nit is given by equation (jH]). In the case of 
a Gaussian distribution of the couplings we find for dt the equation derived in |21j . 
Starting from an initial configuration with do ~ 0, the system is said to be chaotic when 
the Hamming distance satisfies dt > for large times t ^ oo. In the CH-phase two 
paths that are initially close to each other diverge for t — oo. 

We consider a stationary ansatz dt = d and rrit = m in equation (!5T|) . When d > 
we say that the system is chaotic. Because the transitions are continuous we can study 
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the bifurcations around the d = solution. We find the following equation for the 
inverse transition temperature (3* to the CH-phase 



oo oo min(fc°"Sfc'n) 

1= E E^'" E 

fcout>0fcin>0 u=0 

exp {f3*{x + l)) 



dxdy 



exp {p*{x-l)) 



2cosh(/5*(x + 1)) 2cosh(/5*(x - 1)) 

pm 




6{x-2v + t 



-1) • 



(52) 



oo equation ([52]) reduces to T = 4 e^^/^v™ ^ g-2/3>mj gg^^^ ^ 

the different phase transitions are shown. For p large enough, and decreasing the 
temperature starting from a large value we obtain subsequently the following phases: 
P-phase, CH-phase, the chaotic part of the F-phase and the non-chaotic part of the 
F-phase. 



8. Neural Network on a Scale-Free Graph 

The interactions between neurons in organisms are most of the time asymmetric. 
Introducing asymmetric couplings in models for neural networks increases the biological 
realism of the models under study. That is why in [T2| [22] the Hopfield model was 
generalized to include asymmetric couplings. 

We add some more realism to the model by defining the neural network on a graph 
with a given degree distribution. Many real-world networks have a degree distribution 
of the form p{k) = ak~'^ , with a a normalisation constant. These are called scale-free 
graphs. One example is the network of brain activity which has scale- free features |24j . 
In |25] neural networks on scale-free graphs with only symmetric couplings were studied. 

We consider a neural network on a fully asymmetric graph with the following 
distribution of indegrees and outdegrees [26] . 

p{k'\ A;-*) = Aak~^,6 {k'\ fc^^) + (1 - A)a'kr^Xnt ■ (53) 
The correlation factor A in ( l53i) denotes the fraction of sites where the number 
of connections entering and leaving the site are equal. Real-world networks have 
correlations between the indegrees and outdegrees [26]. This correlation between the 
degrees will turn out to have much infiuence on the performance of a scale- free neural 
network. For fixed A we will change the average number of interactions by increasing 
the lower bound b: p(fc™, /i;°"*) = for fc" < b and < b. 

We take the strengths of the interactions Jij according to the Hebb rule: 

/i=i 

with the C,i G {—1,1}, uncorrelated patterns drawn from the probability distribution 
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The network has p patterns = {^l, ^f, ■ ■ ■ , ) on each site. Because the Jij are not 
i.i.d.r.v. variables we can not use the equations f|T71) and f|T8|) . 



8.1. Glauber Dynamics 

We first derive the recursive equations for the marginal distributions when the variables 
evolve through a Glauber dynamics. From these equations we determine the phase 
transition from a P-phase to a retrieval phase (R-phase). In the R-phase the network 
can recover a stored pattern while in the P-phase the noise is too large to retrieve a stored 
pattern from a distorted signal. To calculate the mean of the cavity equations ( fTOj) over 
the quenched variables, it is necessary to define the sublattices = {i G V\$,i = $,}. 

The averaged path probabilities P^(cr(*)) and (c''*^), on the sublattices are defined 



as 



eh 



(56) 



(57) 



When the graph is drawn from an ensemble defined by a degree distribution of the form 
p(/i;™, /c"*^', /c*^^™) = p{k^'^,k"^^)6{k^^"^), such that there are no symmetric couplings, we 
get the following recursive equation for P^ 

Cout [k^^) 



Pi (a-O 



oo 



p{k'' 



Cout 



n EE 



(o))n 



exp /5a(s+l)E 



2P 



0<f'<A:'" 



s>0 



2 cosh 



i3k 

0<f'<fci" p ' 



f'[S 



(5^ 



In the above the Cout(^") is the average number of directed edges leaving a site, given 
the indegree /c"^: Cout(fc™) = J2k°^^ Pi^'^^^W^)^"^^ ■ For the averaged path probability on 
the original graph we obtain analogously 



n 



2P 



(o))n 



exp [/3cr(s + 1) Eo<£'<fc- ^^^'(s 



s>0 



2 cosh 



(59) 

We will use the notation P^(a{t)) = | (l + "^|('f)o"('f)) with superscript a = d or 
a = real. The magnetisations and m^'^'^^ evolve in time according to 



2 Yl 



tanh I P 



Cout 



l + (7,m| ft - 1) 



(60) 
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with the average connectivities 

Cd{k"') = coutik"') = (1 - A)cont + A , (61) 

Creal(/£'°) = Cout • (62) 

We simplify the equations with the condensed ansatz m|(t) = ^^m^(t). This ansatz 
assumes that the spins (= neurons) only have a finite overlap with the first pattern. 
The overlap m^{t) evolves according to 



m 



OO / 7 i 



Cout 



(63) 



with 

r=0 s=0 

"1 + m^(t- 1)" 



1 - m^(t - 1) 



tanh 



P{2s + 2r - fc>) 
P 



(64) 



When calculating numerically the sum in the degrees k^^ in equation fl63|) we have to 
introduce a cutoff K. We will bound by two values, mf < and > m^, with 
and mj^ defined through 

K 



Xt) = Yl p{^'')'^^^^M{¥ 



k"'=b 
K 



Cout 



m 



Cout 



sign (m^(t - i)) Y Pi^' 



Cout 



(65) 



(66) 



In equation (IMl) we used that M(oo) = sign (^m^{t — !))■ Because we have a power-law 
decay of the degree distribution (and not an exponential decay) it is important to take 
the cutoff K into consideration when we want to know the asymptotic behaviour of the 
neural network for K ^ oo. The macroscopic observables will converge much slower to 
the asymptotic value ^^ = 00 when the degree distribution is power law. For finite K 
the time evolution of equation ( l65i) is confirmed by Monte Carlo simulations. 



8.2. The Retrieval State 

When we consider a stationary state of the form, m{t) = mit — 1), we get the following 
equation for the critical inverse temperature 0^ of the P to R transition 



l = ^p(A;-)^i^3.^(/5R }^-^ 



with 



A{r, k 



Cout 



r=0 s=0 ^ 



(67) 



s 
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(2.-nunh('^''^- + ^'-*'»). (68) 

Below the critical temperature Tr = 1/ m^'^^^ > such that the neural network can 
retrieve a stored pattern from a distorted initial configuration. To calculate numerically 
equation (l67j) we take a distribution p(/c™) = ak^ with A;™ G [6, ■ ■ ■ , K] and zero for 
other values of /c"^. We introduce the lower critical value P^, and the upper value 
through 

1 = X:p(n^^(A^,n, 



K 



To derive equation (!69|) we used that C) has the following asymptotic behaviour 
for C — ^ oo 

A{(5,C) ^ J—VC . (70) 

V Tip 

The asymptotic value of the Riemann sum in the second term of fl69p can be calculated 
using a series that converges exponentially in K [2^ . In figure H] we compare how 
and converge to their asymptotic value for K ^ oo. The upper bound clearly 
converges fast to the asymptotic value. The lower bound on the other hand converges 
very slow to its asymptotic value. When we would have bounded the critical value (3^ 
from below with j3^, which we do usually for Poissonian graphs, we would have obtained 
a bad estimate of the critical temperature f5^. Therefore, we estimate in figures [5] and [6] 
the critical temperatures with the upper bounds T^. There we plotted the critical 
temperatures as a function of, respectively, the exponent A and the correlation factor 
A. The retrieval phase increases with the exponent A, which is expected because the 
mean connectivity of the graph also increases with A. Increasing the correlation A 
between the indegrees and the outdegrees on scale-free graphs has a positive effect on 
the performance of the network. The neural network becomes much more tolerant to 
noise and can retrieve considerably more patterns when A increases. 

8.3. The Chaotic Phase 

In this subsection we determine the CH-phase of the neural network. In [22] this was 
done on a Poissonian graph using annealed methods. We consider two systems on the 
same graph undergoing the same thermal noise through the coupled dynamics of 14.21 
We find the following Markovian process for the single time marginals Pt^'- 



H4 



0<f'<fc'" 0<f'<fc'" 



(71) 
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Figure 4. Neural network on a scale-free graphs: The bounds f3^,(3^'^ (upper lines) 
and /3j^,/3f^ (lower lines) on the inverse critical temperatures /3^,/3'^^ as a function of 
the cutoff K. The bounds on (3^ are calculated for the model parameters X = 2, p = 3, 
6 = 4. The bounds on are calculated for A = 2.2, p = 3, b = 4. The upper bounds 
saturate much faster than the lower bounds. 
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X 

Figure 5. Neural network on a scale- free graphs: The critical temperatures (solid 
lines), T'^^ (upper dashed lines) and T™ (lower dashed lines) as a function of the 
exponent A (see equation (|53p ) for a different number of patterns p. The minimal 
indegree is b = A and the correlation factor A = 0. The R-phase is located below the 
solid line. The neural network is chaotic between the dashed lines. 

where Wc is the transition probability defined in ( l28l) . We parametrise the single time 
marginals P^^ with the magnetisations m^ i(t), m^ 2(^) and the Hamming distance d^{t) 
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Figure 6. The critical temperatures of the retrieval state of the neural network as 
a function of the correlation factor A of the scale-free graph for a different number of 
patters p are compared. The ensemble of scale- free graphs has the parameters A = 2.8, 
& = 4. The critical temperatures are estimated with for K = 2000. The R-phase 
increases considerably with the correlation factor A. 



at time step t: 



- [1 + (J m^,i(t) + r m^^2{t) + ar (1 - 2d^{t))] . 



(72) 



The Hamming distance evolves according to 

oo 



n 



1 + aum^^^iit - 1) + Tim^^^2{t - 1) + (1 - 2d^^(t - 1)) cr^r^ 









exp 


5E,(sr) 




2 cosh 








2 cosh 


>E,.(^'-,,) 



(73) 



The magnetisations m^ i(t) and m^2(^) behave according to equation fl60l) . We use the 
condensed initial conditions, m^i(O) = ^^m(0), m^2(0) = .^^m(O) and (i^(0) = (i(0). 
These initial conditions imply that the initial configurations have a finite overlap with 
only the first pattern. From the time evolution (1711) we find that for the condensed 
initial conditions m^^i(t) = C,^m{t), m^^2{t) = $,^m{t) and d^{t) = d{t). The evolution of 
the overlap m{t) is given by equation ( l63l) . For d{t) we get 

oo fc'" 



d{t) = ^ p{kn E 




(rf(t-l)) 



c/(t- 1))" 
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(fcm_„)(p_l)„(p_l) 

Rl=0 _R2=0 



77, 



c/xfiy exp 



-n{p- 1) + 2i?2 + x) 



exp 


} i\y\ 


+ 2i?i - (A;i°-n)(p- 1)) 




2 cosh 






-A:i-(p-l) + 2i2i + 2/22) 



exp 


- (- 


y\ -2Ri + (A;'° 1)) 




2 cosh 




+ (A;i>^ - 2n)(p - 1) - 2Ri + 2i?2)] 



v=0 



n 

V 



1 +m(t- 1)/(1 i))y 



l-m(t-l)/(l-d(t-l)) 



5 [x — 2v + n) 



w)=0 



t'^ -n 
w 



""5 (y - 2m; + A;"' - n) . 



(74) 



In equation ( 17^ we summed subsequently over the following variables: the indegrees 
k^^, the number of neighbouring spins n with cr 7^ r, the number of neighbouring spins 
V with 0" = r and a = 1 and the number of neighbouring spins w with a ^ t and 
(T = 1. The summation variables Ri and i?2 are the number of non-condensed patterns 
on the neighbouring spins with, respectively, a ^ t and cr = r that are equal to the 
corresponding pattern on the original site. The complex function /(i?; x) used in 
equation fITil) equals 

r.27r 



1 
2^ 



dio exp [zu;i?] 



cos 



a; 



exp 



—ix arctan 



sma; 

1 + cos LO 



(75) 



It is possible to find an equation for the transition temperature (5'^'^ to a CH-phase with 
the stationary value d > of (i(t) by expanding the left side of ( 17^ around dit — 1) = 0: 



1 = ^p(yti-)yti-^:^ — Le{^, t'^^m) 



with 



dx 



v=0 



Cout 

A;"^ - 1 

V 



1 + m 



fc'°-l-r> 



5 (x - 2t; + A;"^ - l) 



p-l (fci"-l)(p-l) 

E E /(^i>p-i)/(^2,(p-i)(fc''^-i)) 

i?i=0 ^2=0 
'/5 



exp 



(-(A;'°-l)(p-l) + 2i?2 + a;) 
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exp 


^ (1 + 2i?i - 


{P- 






2 cosh 




f 1) - 


1)- 


f 2i? 


1 + 2R2) 



exp 


f (-l-2/?i + (p-l))' 




2 cosh 


^ ((a; - 1) + (-A;- + 2)(p - 1) - 2i?i + 2i?2)] 



The asymptotic behaviour of B for C ^ cxd is given by: 
i3(/5,C,m) ~ exp(C^(m)) , 



(77) 



(7^ 



with 

\l/(m) = log 
and 



'1 — m) + (1 + m) exp 



{P - 1) log 



1 + exp 



2x 



X , (79) 



X = - log 



2m 
P 



m 



—m + 



2m 



j — (m — 1)(1 + m) 



1 + m 



(80) 



When |m| > we have that \E' < 0, hence, the series converges exponentially for large 
C. When m = we have that \& = 0. The asymptotic behaviour of B is then: 

S(/?,C,m)~^('./A^2-P+iV| - ^ )|l + 2r-(p-l)| = ^. (81) 



C 



We define the upper bounds and Pi on the inverse critical temperature /? to the 
CH-phase for m = as: 



K 

E 

K 



Cout 



1 = ^ p(k^)t-^^B{(3^^, t\ o)+c Yl p(kn^ 



Cout 



c(A;'°; 

Cout 



(82) 



(83) 



The convergence of (5'^ and ^if" to their asymptotic value fi'^'^ in function of K is plotted 
in figure m Because (5'^ saturates much faster we used this value in figure [5] to estimate 
The R-phase contains a part with d > 0. The bounds /5™ and /5™ on the inverse 
critical temperature /3™ of the transition from the chaotic part of the R-phase to the 
non chaotic part of the R-phase are calculated with substitution of respectively m\ and 
into equation (175]) . The value of /9™ is a good approximation to In figure O 
the complete phase diagram of the neural network with the P-phase, the R-phase and 
the CH-phase is presented. The chaotic region of the neural network is enclosed by 
the dashed lines. This region is larger for odd values of p. The R-phase and CH-phase 
become smaller when p increases and the non chaotic part of the R-phase dissapears 
when 7 increases. The CH-phase is larger for odd values of p than for even values of p. 
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9. The Stationary Solution 

In this section we develop an algorithm to calculate the marginals of the stationary 
solution. This algorithm would be the equivalent of the belief propagation equations 
( ]B.14[) . Equations ( fTOl) constitute an algorithm with a linear computational complexity 
0{N). But because the algorithm scales as (9(2*) in time t we can only follow the 
dynamics for a small number of time steps. We are interested in the marginals of the 
stationary state when t oo. To solve this we introduce some assumptions on the 
distribution P(cr° '*). 

9.1. The One-time Approximation 

The first simple approximation one can make is to neglect the correlations in time: 

In general (see for example equations (l37j) . (l39l) and (HTl) ). when solving the dynamics, 
we have equations of the following type 

P{a'-'\e'-') = T [o'-'r-'- {P,, J,}.=i...) • (85) 

We close these equations with (l84l) using 

P\o{t)\e)= \mi y P{a'-'\e'+^-'), (86) 



s—*~oo 

„s..t~l 



with 6'*"* the constant vector with components 6. When we insert equation (!84l) in the 
left hand side of fl86l) we find 



P*{a{t)\e)= hm 5^ . (87) 



S^ — OO 

„s..t-l 



This approximation leads to the correct stationary solution in the case of models defined 
on fully symmetric or fully asymmetric graphs. This makes us curious to see how 
good this approximation would work for models defined on partially asymmetric graphs. 
Because we neglected correlations in time we can not expect a good description of the 
spin glass phase. 

9.2. Density Evolution Equations in the One-time Approximation 

We apply the approximation given by flMl) and fISB]) to the equations (jSTj), fl5Ul) and fHT]) : 

(88) 



t 



L 2 cosh [/5m<^ 

[exp [Pu (t( 
^ 2cosh[/5n] 
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We insert ( 155]) . and (pi in the update functions T given by equations 

f HOj) and fH2l) and close the equations through fl86l) . We define the function u 
Wfcin fcsym Jf} ; ^) tfirough the explicit solution u of the implicit equation 



u 



in 



exp 




2 cosh 


/5 (EtV^^^ + ^) 



exp i/3 xi?rr ^^^^ ('^^^) + EEfc^ym+i ^^^^ 



n^lT cosh J^r)] 



exp [/5ru] 



(91) 



In the sequel we use bimodal distributions of the form fj43|) . We then have two type of 
messages propagating along symmetric links: u[ff) = u{Joa) = u"^ . Using the one-time 
approximation in (1371) we get for the density of the fields propagating along directed 
edges 



oo min(fc°"',fc'") 



Cout Csym 



5 [n - Wfcin^fcsym ({n^/, J^/}^,^-^ ; 0)] . (92) 

The density of the fields propagating along the symmetric edges follows from the 
equations (l39l) 



fcin>0 fcsym=0 



-sym 



fj f dJ,R{Ji) I du,W\u,) 



JJ I dJiR{Ji) I dujdujW'y"'{ui,u} 

£=1 



(93) 



b\u — Wfcin_l ;jsym_l }^' = 1. 1 ? ""'-)] 

5 [n"*" — Wfcin_l fcsym_l ({m^', J^'}^/^]^ j,in_l ; l) ] • 

Substitution of the one-time approximation (|90|) in the equation for the distribution 
of the real marginals (HTj) gives 

oo fc'n fc=y™ „ „ 



fcin>0 M=0 



JJ /" dJeR{Je) ! due (n^) 5 [n - Wfein^^syn. ({u^/, J^/}^,^^ ; 0)] (94) 

Because of the analogy with the equations (IB.12P we call the equations (|92|) and (|93|1 the 
density evolution equations in the one-time approximation. We see that now we have 
two densities instead of one: the density for the fields propagating along symmetric 
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edges and the density for the fields propagating along directed edges. The equation for 
the marginals on the original graph is the equivalent of flB.lSP for the symmetric 



model. Another important difference is that the update function W^in ;,sym, for partially 
asymmetric graphs, is not explicitly known. Instead we have to solve the implicit 
equation (19T1) . The equations ( l92l) . ( l93l) and ( IMl) constitute an algorithm that generalizes 
belief propagation to graphs with asymmetric bonds. 

9.3. Fully Symmetric and Fully Asymmetric Cases 

For graphs with exclusively symmetric or asymmetric links, we determine the explicit 
form of Z///;in fcsym . For a fully symmetric graph, the equation ( 191]) admits the solution 

Wfcin_i_fcsym_i JJ^^^ j.sym_i ] 9) = 6 + (3~^ ^ atauh [tauh {(3Ji) tanh {(3ui)] , (95) 

where we used ue(JeT) = .J^t + u^. For fully asymmetric models we have the solution 

W,. ,.y. ({u., J,},^,.,,. ; 0) = }^ a log O , , r.,v, 



2/? ^ ^ I ^ 2 cosh [/5 (E, -hr,)] '""^ 



13 ^ nui, 
1=1 



(96) 

From (!95|) and (!96|) it follows that the one-time approximation gives the correct results 
for Glauber dynamics of the Ising model on fully symmetric and asymmetric graphs. 
Now we check how good the equations fl92l) . fl93l) and flMl) are for models on partially 
asymmetric graphs. 

9.4. Results and Comparison with Simulations 

We numerically solve the equations (1921) . (!93|) and (!94l) with a Monte Carlo integration. 
The unknown distributions W^, W^^'^ and 1^''''^^ are represented as populations. The 
procedure is also known as population dynamics [7] . In figure [7] the magnetisation is 
plotted as a function of the temperature for an Ising model without bond disorder on a 
Bethe lattice. The degree distribution is then: 

p(A;'^ Fy'") = 5{k''^'^ - C)6{k"' - D)p{k°''') . (97) 

Because there is no disorder the distributions W^lu) and iy^^™(n+,n~) are delta 
functions. Figure [7] tells us that the theory and the simulations are in good agreement. 
For C = 0, 3 the results coincide while for C = 1,2 there is a small deviation. In figure 
|8]we checked how good the approximation is on graphs with fluctuating connectivities. 
We simulated an Ising model without bond disorder on Poissonian graphs drawn from 



the ensemble defined in [Appendix A[ Despite the fiuctuations in the degrees we find a 
good agreement between theory and simulations. This confirms that for models without 
bond disorder, the one-time approximation works very well. In figure [9] we plot the 
magnetisation on a Bethe lattice with bond disorder. We see that the difference between 
theory and simulation increases with the bias p in the bonds. 
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Figure 7. The magnetisation m as a function of the rescaled temperature T/{DJq) 
for a ferromagnet (i.e. p = 1) defined on a Bethe lattice. The indegree D of the graph 
equals 3 and each site is incident to C symmetric bonds. The simulations (markers) 
are mean values of 20 runs on graphs of sizes O(IO^). The theory (lines) follow from 
solving recursively the density evolution equations in the one-time approximation. 




ao.4- 



T/(cJ„) 



Figure 8. The magnetisation to as a function of the temperature T/[cJo) for an Ising 
model without bond disorder on a Poissonian graph with mean connectivity c — 3 and 
different fractions of symmetric edges e. The lines are obtained by population dynamics 
from the density evolution equations in the one-time approximation for populations 
of sizes O(IO^). The markers are the average results from 20 runs with the heat-bath 
algorithm on a graph instance of size 0(10^). 
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Figure 9. The magnetisation m as a function of the bias p in the couphngs (see 
equation (|43|) ) at a temperature T/{DJo) = 0.5. The Ising model is defined on a 
Bethe lattice with the degree distribution defined in (|97)) . The indegree D equals 3 
and results arc shown for various values of the symmetric degree C. The simulations 
(markers) arc the mean values of 20 runs on graphs of sizes O(IO^). The theory (lines) 
arc results from the density evolution equations in the one-time approximation using 
a Monte Carlo calculation with populations of 0(10^) fields. 



9.5. Bifurcation Analysis 

We determine the P to F and the P to SG transition lines for Ising models with bimodal 
distributions using a bifurcation analysis around the paramagnetic solution. First we 
note that the equations fl92l) and fl93l) admit the solution: 

W'^iu) = 6{u) , 



W'{u+, u~) = dA W^{A)6{u+ - A)6{u~ + A) . 



(98) 
(99) 



Indeed, when we insert ( 198|) and (!99|) in (!93|) we get for W^{A): 



w-(A) =Y.Y. 

fcin>0 fcsym = 



^sym 



n 



dAfW^{At 



0<f<fc=>'™-l 

5 [a — ^fcsym_l fcm_l ( { j^/ = 1 . .fesym _ 1 ) ] 

with ^fcsym_i ;-in_i thc expllclt solution of 



(100) 



--±1 



exp 




cosh 


PJo 





cosh 



(101) 
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The solution (1981) and fl99p represents the paramagnetic phase. For a symmetric lattice 
with /c'° = k^^"^ we have Ak<>ym_i ^in^i = Jq. For k^^"^ = we obtain 



1 I exp 



2P 



=±1 



T cosh 



(102) 



For partially asymmetric graphs W^{A) results in a more complicated form. 

To calculate the P to F and P to SG transitions, we expand the equations ( !92l) . ( l93l) 
and ( IMI) around the solution given by ( l98i) and ( l99l) . We use an expansion in a <^ 1 
with: 



and 



(a) = / du+du~W{u+, u~) {u{a) - aA)" ~ 0(a") 



rfM+c/M-l^(M+, M-) (m(+) - A)" (m(-) + A)" ~ C(a"+™) . 



(103) 



(104) 
(105) 



Substitution of ffT03|) . ffTOij) and ffTOSD in ([92]), ([93D and ([941) gives, up to linear order in 
a, the equation 



m 



pM 



m 



m 



with 



M 



Mdd Mds 



(106) 



(107) 



The critical temperature Tp of the P to F phase transition is given by the equation 
\\\~^(Tp) = p with A(T) the eigenvalue of M with the largest modulus. The elements 
of the M matrix are 



fc°"'>Ofcin>0 fc=>""=0 
;.sym 



Cout '^sym 



(108) 



ds 



oo oo mm(A;°"Sfc'") 
A:°"'>Ofcin>0 /c=ym=0 



Cout C-syjn 

(ri)d-tanh {(3Ai) (r)/ 



y^sym 



J] dA,W{Ai 



g |; p(fcin,fcsym)fcsym ^^^^^^^ 



(109) 



--sym 
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Figure 10. The phase transition hnes between the P-phase. F-phasc and SG-phase on 
a Bethe lattice with fc'" = 3 and k^y"^ = C. The hnes are calculated from the equations 
for the critical temperatures derived from a bifurcation analysis. The markers are 
computed with the population dynamics method on the density evolution equations in 
the one-time approximation. 



II dA,W{A, 



1 - (r). 



fcin>0 fcsym = o 



--sym 



■J 1 



{Tk^y n-1 ) , 



The averages {■)d and {■)s appearing in the matrix elements are given by 

(/(r,r, A))rf = ^Wfcsym^fcin (l,T,T;0,A(^{Ae}^ f^syn.^^ f{T,T, A) , 
T,T 

and 

(/(T,r, A))s = ^ o-Wfcsym_i,fcin_i (|cr,r,T; Jq, ^ ({^Ji,fcsym_i)) /(T,r, A) 



(110) 



:iiii 



(112) 



r,T,(T 



(113) 



in the weight Wfcsym ^in. The weight Wfcsym is expressed as 

Jo 



exp 



>VfcBym_fcin(a,r,r;6', A) = 



2 cosh 



/3 JoEf=i^«+e 



exp 



/5^(E'17^^^^ + sign(e)^ 



oxp 



/3a(joEf=Vf+e) 



2 cosh 



[/3(JoE£ 



exp 



(5t ( E£T ^^^^ + sign {B) A 



(114) 

For symmetric graphs the bifurcation condition |A|~^ = p of fll06p simplifies to 
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ml = pMssTnl, which leads to the condition 

.tanh(/3^Jo) ^ pik^n = 1 ■ (115) 



Equation flllSp gives the P to F transition line on a symmetric graph [28]. For fully 
asymmetric graphs we get mf = pM^^m'l: 



(116) 



This is the same condition we found in (l49l) . 

To determine the P to SG transition we expand up to order with J duW'^{u)u = 
and J du W^{u'^ , u~){u{a) — a A) = 0. We then arrive at 



TTln 



mo 



Q 



rrin 



, Q 



Qdd Qds 
Qsd Q ss 



(117) 



with 

Qdd = 



oo oo 



mm(fc°''S/fc'") 



Y[ dA,W{A, 
1=1 

oo oo min(fc°"Sfc'") 

E E E 



2 



1 - (r). J 



:il8) 



Cout Cc 



sym 



(ri)d-tanh (/3Ai) (r)d 



(119) 



oo fc'" 



<3»=E E 



p{k^,k''^"')k 



sym 



,sym 



^sym 



(ri),-tanh (/3Ai) (r), 



oo fc'" 



^(^^m^^sym^^sym 



^sym 



Q«.= E E 

Again we find the correct bifurcation condition 

tanh^ iPscJo) E Mfe^^-) ^^3,^ ^ 



n 2 



(120) 



:i2ii 



fcin>0 



(122) 



for symmetric lattices [28] . In figure [TO] the phase transition lines between the P-phase, 
F-phase and SG-phase are shown for a bond-disordered Ising model with a bimodal 
distribution. The P to F and P to SG lines are obtained using the bifurcation analysis 
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Figure 11. Critical temperatures computed from a bifurcation analysis (lines) within 
the one-time approximation are compared with simulations (markers) as a function of 
the fraction of symmetric edges e. The results shown for the P to F transition are for 
an Ising model on a Poissonian graph without bond disorder. The results shown for 
the P to SG transition arc for an Ising model on a Poissonian graph for p = 0. 



while the F to SG transition is found through population dynamics. In figure [H] we 
present the P to F transition line for a ferromagnet (p = 1) and the P to SG transition 
line obtained through the bifurcation analysis on a Poissonian graph as a function of the 
fraction of symmetric edges e. The markers from simulations are in agreement with the 
theory. The simulations are performed through the method of Binder cumulants, see 



[19] . In Appendix C we give some details on how we derived these critical temperatures. 
We see that the SG-phase dissappears when the asymmetry in the graph increases. 



10. Conclusion 



In this paper we applied the cavity method to study the dynamics of spin models on 
a given graph instance. We derived a set of effective equations which describe the 
dynamics. Solving these recursive equations can be seen as the equivalent of the belief 
propagation algorithm known from inference problems or decoding algorithms. Just 
like the latter, we expect these equations to be exact on a tree. The main difference 
with statics is that path probabilities, instead of stationary probabilities of single spins, 
are propagated along the edges of the graph. We took the average over an ensemble 
of graphs to find the recursive equations describing the dynamics of Ising models on 
typical graphs drawn from this ensemble. These equations generalize the result of 
[T5] to graphs with arbitrary degree distributions. The macroscopic evolution of the 
system is given as a function of three mean values of path probabilities: the one of the 
probabilities propagating along directed edges, the one of the probabilities propagating 
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along symmetric edges and the one of the marginal path probabilities of the spins on 
the original graph. On a Poissonian graph the three recursive equations for these path 
probabilities reduce to one and we find back the result derived in [15]. We solved the 
following problems on the basis of these equations. 

The evolution of macroscopic observables like the magnetisation is calculated for 
a small number of time steps because the computational complexity of the derived 
equations scales exponentially in time. These results are compared with simulations 
and both methods are in good agreement. This confirms that the derived equations are 
exact on graphs with a local tree structure. 

The effective equations simplify on fully asymmetric graphs. The effective dynamics 
becomes Markovian and the stationary distribution is derived. We studied the phase 
diagrams of an Ising model on a Bethe lattice showing the existence of a paramagnetic, 
ferromagnetic and a chaotic phase. We also examined a neural network with Hebbian 
couplings on a scale-free graph where we found that correlations between indegrees and 
outdegrees increase the retrieval properties of the neural network tremendously. 

We introduced an heuristic method to derive the stationary distribution of Ising 
like models on partially asymmetric graphs by assuming that the path probabilities 
factorize in time. The equations are closed such that we recover the correct stationary 
distribution for models on fully symmetric and fully asymmetric graphs. This method 
is approximative. We found a set of equations that generalize belief propagation to an 
algorithm that computes the marginals of the stationary distribution of Ising models 
on graphs with asymmetric couplings. The evolution of macroscopic observables and 
the phase transitions are obtained with this method and compared with simulations. 
Without bond disorder the theory and the simulations give almost identical results. 
When introducing bond disorder the theoretical predictions start to fail and a more 
refined approximation is necessary in that case. 

In future work we will try to get better algorithms that derive the marginals of the 
stationary solution by using a systematic method that minimizes a functional like the 
KuUback-Leibler distance. Work along this way is in progress. Furthermore, it would 
be interesting to analyze the dynamics of a neural network, for example given in |29j . 
with the cavity equations. 

Appendix A. The Poissonian Ensemble 

The Poissonian ensemble Gp{c, e) is defined through the following probability function 
Pp(C; c, e) = n (^^(^^^' 1) + (l - ^) 0)) 

i<j 

n {e6{c,,- c,,) + (1 - e) (^5(c,,; 1) + (l - ^) 0))) . (A.l) 

i>j 

The mean connectivity is given by c while e is the fraction of symmetric edges. This 
is the same ensemble as used in [20]. In the recursive equations f|T7j) . f|T8|) and f|T9l) for 



We write the Kronecker 6 as 



AT 



The Cavity Approach to Parallel Dynamics of Ising Spins on a Graph 34 

dynamics we need to know the degree distribution. We calculate the degree distribution 
Pp, defined through 



(A.2) 

(A.3) 
(A.4) 
(A.5) 



j 

"'0 



2- Ay 

271 



exp 



j 

iuj'C^^ Cji — k 



Substitution of the above integrals in (lA.2p leads to 

doj duj' d( 



e-e"+'^'+'^ + e ( 1 - 
N 



2n 2n 2tx 
c 



N 



exp [-lujk^'' - iu'k'"'' - iCk""^""] 

+ e'^') + (1 - e) ( 1 



l-^)-fl-- 

'N \ N 



c x2 
N 



We keep only the O (A") terms: 

p,{k-,k^^\k^n fd^d^'d^ 



exp 



ce e 



27r 2n 2tx 



exp {-iujk''' - iuj'k°''' - iCk'^""] 

(1 - e) c f + e"'") - 2 (1 - e) c 



(A.6) 



When we expand the exponentials in flA.6P we get 

du doj' d( 



Pp{k''\k°''\k'^'^) 



2tt 2tx 2tx 
J exp [mi {uo + uo' + Q] 



exp [-lujk"' - tuk"""' - tCk'''"' - ce - 2 (1 - e) c] 



711 



rill 



"2 



»^2 



exp [in2Uj\ 



^13 



exp [zriso; 



Integrating over gives us 



27r 27 



exp [-iujk''' - tu'k"''^ - 2 (1 - e) c] 



Y: (c(1 - e)r [^("^ + ^""^)"] 5^ (c(l - e)) 



exp [iim + k'y'^)uj'] 



n2\ " ns! 

n2 na 

with Pp(/c; c) = exp(— c)c'^/fc!. Finally we integrate over uj and u': 

Pp{t\ A;-*, A;^^-) = pp (fc^^^; ce) Pp (A;'-^ - A;^^'"; c(l - e)) Pp (A;-* - A;^^-; c(l - e)) 

When we sum out the out degrees we finally obtain 

g-c^fc'" / ^in 



^inl 1 ^sym 



(A.7) 



(A.8) 



The Cavity Approach to Parallel Dynamics of Ising Spins on a Graph 35 

For the Poissonian ensembles the equations for the dynamics ( I17p , (1181) and fll9p simplify 
because of the property 

Lsym 

^in ^sym Sym ^sym 

(A.9) 

Indeed, repeatedly applying dAM on ( ITTl) . ( fT8|) and ( fT9l) gives P^^'^^ = P^^^ and when 
Qi..t _ Qi..t. p _ psym ^j^^ three equations reduce to one equation. 



Appendix B. Stationary State for Glauber Dynamics 



By using the tools of equilibrium statistical mechanics we determine the stationary 
behaviour of a model evolving through Glauber dynamics determined by the equations 
([T]) and ( !22l) . We consider the case where the field equals 

hi{s) = hi{(T{s - 1)) = ^ JijCijaj{s - 1) , (B.l) 

3 

with Jij = Jji and Cij = Cji. The stationary state is given by Pst(c) ~ exp {—PH{cr)), 
with H{(t) equal to 

^('^) = E (2 cosh . (B.2) 



1=1 



This system has the same thermodynamic behaviour as a model with the two-spin 
Hamiltonian 



(B.3) 



h3 



Performing the standard replica method or the cavity method (see for example [30] ) 
we get the following set of selfconsistent equations for the distribution W of the cavity 
fields M, V and w, within the replica symmetric approximation: 

W'{u,v,w) = Y^^ l[{dJiR{Ji)) / l[{dudvidweW'{u,,v,,w,)) 

J n ^ ^ D 1 n 1 



fc=0 



fe-i 



u 



w 



ui, ve, We 
ue, ve, We) 



k~l 



i, Ui, Ve, We 



with 
U{J, u, V, w 

V( J, u, V, w 



) = — Y^ ^ log I 6xp ( Jst + Jst + us + vt + wsi)] 
— t log < exp ( Jst + Jst + us + vi + wsi)~\ 



(B.4) 

(B.5) 
(B.6) 
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>V( J, u, V, w) 
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^.i [ ~s,i J 



36 



(B.7) 



The distribution determines the density of the fields (m^^ u^P'') along the edges 
of the graph. These fields correspond with the marginal probability distribution pj^''^ on 
the cavity graph G^^^: 



exp 






exp 


(j) (j) (j) 
ul a + VI T + wl ar 



(B.8) 



On a given graph instance it is possible to calculate the fields {u\''\vl''\ w^^) through 
the belief propagation algorithm. In the n-th time step of the algorithm, messages 
(Uilij,vl'^j,wl'^j) are propagated along the edges of the graph. The update equations 



of these messages are given by 

E7y( T ^,>('") ^ 



(n+l) 



fn+l) 



■ . _ \^ ■])( T 1,/") ^ 



w 



(B.9) 
(B.IO) 
(B.ll) 



The fields {uf\v\^\w\^'^) are given by the solutions of (iRQl) . (iRTnjl and (IRTTI) for 
n ^ oo. We call the equations (lB.9p . fIB.lOp and (1B.11|) the cavity or belief propagation 

can be found in equation 



7!- ■ Ul - 
1 t^p ^i^jl 



equations. The density of the messages (Uj-J^ 
(IB. 41) . That is why ( IB. 41) is also called the density evolution equation. 
When we choose W{u, v, w) = 6{w)6{u — v)W'^{u) we find 

" p{k)k "^^^ "'^"^ 



W%u) = / \{{dm{J,)) j \{{du,W\u,)) 

k=0 £=1 l=\ 

fe-1 

— atanh (tanh(/3J£) tanh(/3M^)) 



(B.12) 



The distribution of the real fields u of the single site marginals X]cr\o- -Rit(''') 
exp [wjO"] are given by the selfconsistent equation 

W^{u) = J2^^ / l[{dJ,R{J,)) / HidueW'iue)) 

1 n ^ J f) 1 J D 1 



fc=0 



e=i 



u 



(3 ^ atanh (tanh(/3J^) tanh(/5M£)) 



(B.13) 



The solution W{u^v^w) = 6{w)6{u — v)W'^{u) corresponds with the solution w\ 



(n+l) 







and u. 



(n+l) 



of the equations flB.9p . fIB.lOp and fIB.lip . The cavity equations 
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reduce in this case to 

mJ^+.i) = atanh (tanh {f3Ja) tanh • (B.14) 

The single site marginals -Pi(cr) can be calculated through 

Ui = lim atanh ^tanh (/^J^) tanh (^Pu'l^^^ . (B.15) 



Appendix C. Simulations 

To determine the critical temperature from the P-phase to the F-phase of Ising models 
on a graph of size we calculate the Binder cumulant B]^ [19]: 

Bn{T) ^ 1 - . (C.l) 

The brackets denote an average over the stationary distribution. The quantity m is the 
magnetisation (Xi/N. In the F-phase 5 = 2/3 while in the P-phase B = 0. The point 
Tp where the Bn(T) lines for different values of cross is the critical temperature. In 
figure [CTI we present I^at as a function of T for a Poissonian graph with e = 1 and c = 3. 
The value Tp for finite system sizes is found to be higher then the theoretical result 
Tp = (c atanh(l/c))~^ valid for 00. 

The critical temperature for the P-SG transition is calculated through An [3T| : 



MT)^\{3-^) (C.2) 

The bar denotes the average over the quenched variables. Above the critical temperature 
TsG the lines for different system sizes join into one line. 
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